Strain Hardening of Polymer Glasses: Entanglements, Energetics, and Plasticity 
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Simulations are used to examine the microscopic origins of strain hardening in polymer glasses. 
While stress-strain curves for a wide range of temperature can be fit to the functional form predicted 
by entropic network models, many other results are fundamentally inconsistent with the physical 
picture underlying these models. Stresses are too large to be entropic and have the wrong trend with 
temperature. The most dramatic hardening at large strains reflects increases in energy as chains 
are pulled taut between entanglements rather than a change in entropy. A weak entropic stress is 
only observed in shape recovery of deformed samples when heated above the glass transition. While 
short chains do not form an entangled network, they exhibit partial shape recovery, orientation, and 
strain hardening. Stresses for all chain lengths collapse when plotted against a microscopic measure 
of chain stretching rather than the macroscopic stretch. The thermal contribution to the stress is 
directly proportional to the rate of plasticity as measured by breaking and reforming of interchain 
bonds. These observations suggest that the correct microscopic theory of strain hardening should 
be based on glassy state physics rather than rubber elasticity. 

PACS numbers; 61.41. +e,81.05.Kf,81.40.Jj,81.40.Lm 
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I. INTRODUCTION 

Glass forming polymers are of great industrial im- 
portance and scientific interest because of their unique 
mechanical properties, which arise from the connectiv- 
ity and random-walk-like structure of the constituent 
chains. At large strains, the stress increases as the chain 
molecules orient, in a process known as strain hardening. 
Strain hardening suppresses strain localization (crazing, 
necking, shear banding) and is critical in determining ma- 
terial properties such as toughness and wear resistance. 

Traditional theories of glassy strain hardening P, Q as- 
sume that polymer glasses behave like crosslinked rubber, 
with the number of monomers between crosslinks equal 
to the entanglement length N^- The increase in the stress 
a due to deformation by a stretch tensor A is attributed 
to the decrease in entropy as polymers stretch affinely 
between entanglements. Beyond the initial plastic flow 
regime, entropic network models predict 0: 



(1) 



where Tfiow is the plastic flow stress, is the strain 
hardening modulus, is the inverse Langevin function, 
g{X) describes the entropy reduction for ideal Gaussian 
chains, and L^^{h)/3h corrects for the finite length of 
segments between entanglements. The value of enters 
in h, which is the ratio of the Euclidean distance between 
entanglements to the contour length. 

Entropic network models have had much success in 
fitting experimental data [3, |3, 0, @. However, serious 
discrepancies between the models and experiments are 
revealed in trends with temperature and the values of fit 
parameters @. Entropic models predict Gr = peksT, 
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where is the entanglement density. This is about 100 
times smaller than values measured near the glass tran- 
sition temperature Tg [2j. Moreover, Gr increases mono- 
tonically as T decreases 0, 0] , while any entropic stress 
must drop to zero as T — > 0. Parameters such as Gr and 
A^e, which entropic models assume to be material con- 
stants, must be adjusted significantly to fit data for dif- 
ferent strain states (i.e. uniaxial or plane strain) Q. Fits 
to the shape of strain hardening curves also yield smaller 
values of A'^e than those obtained from the plateau moduli 
of melts ,2J • 

A more fundamental conceptual difficulty with en- 
tropic models is that, unhke rubber, glasses are not er- 
godic. For T < Tg, thermal activation is not sufficient to 
allow chains to sample conformations freely. Rearrange- 
ments occur mainly under active deformation [ill 
and at a frequency that scales with the strain rate p^ . 
In such far from equilibrium situations, the relevance of 
entropy is unclear. In addition, experiments [isl . [3, 
and simulations , [l^ show that the internal energy 

contributes to strain hardening, but this is not included 
in entropic network models. 

Molecular simulations allow a full analysis of the mech- 
anisms of large strain deformation in glassy polymers 
la HH, m m [13, m, m. in recent papers 
have examined the origins of strain harden- 
ing using a coarse-grained bead-spring model [13] ■ As 
in experiments, numerical values of Gr are much larger 
than predicted by entropic models and show the opposite 
trend with temperature @ . A direct correlation between 
Tfiow and Gr was discovered that allowed curves for dif- 
ferent interactions, strain rates and temperatures to be 
collapsed Q. Substantial strain hardening was observed 
for chains that are shorter than and thus do not form 
a network [13, H^- For chains of all lengths, strain hard- 
ening was directly related to strain-induced orientation 
and the rate of plastic rearrangement [17 1. 
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This paper extends our simulation studies in several 
directions. Uniaxial and plane strain compression are 
examined for a wide range of iVg, T and chain lengths. 
Stress curves for all entangled systems can be fit to Eq. 
[Tl The fits confirm the connection between Tfiow and Gr 
[8|, which both drop linearly to zero as T rises to Tg. This 
observation motivates a modification of Eq. [T] Using Tg 
and the fit to a stress-strain curve at one temperature, 
the model predicts strain hardening curves for all T < 
Tg remarkably well. However, as in experiments it 
is necessary to vary parameters such as Gr and N^. in 
unphysical ways in order to fit curves for different strain 
states. 

Direct examination of entropic and energetic contribu- 
tions to strain hardening reveals qualitative failures of 
network models. The rapid hardening at large strains 
that is fit by the Langevin correction in Eq. [T] does not 
reflect a reduction in entropy. Instead there is a rapid rise 
in energetic stress as chains are pulled taut between en- 
tanglements. Variation in the energetic contribution for 
different strain states leads to changes in fit values of iVg. 
For all chain lengths and N^, we find that the thermal 
part of the stress correlates directly with breaking and 
reformation of van der Waals bonds during deformation. 
This provides an explanation for the correlation between 
Gr and Tfiow 

Remarkable shape recovery is observed in experiments 
when highly deformed samples are unloaded and heated 
slightly above Tg [28|. Network models assume this re- 
covery is driven by a "back stress" related to the entropy 
of the entanglement network, and shape recovery is often 
seen as providing evidence for entropic strain hardening. 
Our simulations also show dramatic shape recovery that 
is driven by orientational entropy. However, the magni- 
tude of the associated stress is only of order peksT and 
thus much too small to be significant in strain hardening. 

Changes in microscopic chain conformations are also 
explored. While Eq. [T] assumes affine deformation of seg- 
ments of length Ne, the observed deformation becomes 
increasingly subaffine as strain increases. This reflects 
a straightening of segments between entanglements that 
disturbs the local structure of the glass and increases 
the internal energy. Although uncntangled chains do not 
form a network, they still undergo significant orientation 
during strain @, l25|. For all chain lengths, the thermal 
contribution to the stress is directly related to the orien- 
tation of chains on the end-to-end scale rather than the 
macroscopic stretch [13, . 

The following section describes the potentials, geom- 
etry and strain protocol used in our simulations. Next, 
fits to entropic network models are examined, and an ex- 
tension that incorporates the correlation between Tfio^ 
and Gr is presented. This is followed by a discussion 
of the energetic and entropic contributions to the stress 
and the role of entropic back stresses in shape recov- 
ery. Sections IIII Dl and IHI El explore the effect of chain 
length and orientation and demonstrate the connection 
between plastic deformation and the thermal component 



of the stress. The final section presents a summary and 
conclusions. 



II. POLYMER MODEL AND METHODS 

We emp loy a coarse-grained bead-spring polymer 
model [23] that incorporates key physical features of 
linear homopolymers such as covalent backbone bonds, 
excluded-volume and adhesive interactions, chain stiff- 
ness, and the topological restriction that chains may not 
cross. Each linear chain contains N spherical monomers 
of mass m. All monomers interact via the truncated and 
shifted Lennard-Jones potential: 
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(2) 

where Vc is the potential cutoff radius and ULjir) = 
for r > Tc- We express all quantities in terms of the 
molecular diameter a, energy scale uq, and characteristic 
time tlj = yjma^ /uq. 

Covalent bonds between adjacent monomers on a chain 
are modeled using the finitely extensible nonlinear elastic 
(FENE) potential 
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with the canonical parameter choices 27] i?o ~ l-5a and 
k = 30wo/a^. The equilibrium bond length Iq 2± 0.96a. 
The FENE potential does not allow chain scission, but 
the maximum tensions on covalent bonds for the systems 
studied here are well below the critical value for scission 
in breakable-bond models [3(]| . 

As a means of varying TVe, we introduce chain stiffness 
using the bending potential 



Ubend[0) = fc6e„d(l " COsO) , 



(4) 



where 9 is the angle betweeen consecutive covalent bond 
vectors along a chain. Increasing fcfcend increases the root- 
mean-squared (rms) end-to-end distance of chains Ree- 



The chain stiffness constant Co 



for well-equilibrated jSlJ melt states rises from 1.8 to 3.34 
as kbend is increased from to 2.0uo. The value of 
decreases from about 70 to about 20 over the same inter- 
val [3^ . The key parameter in entropic network models 
is the number of statistical segments per entanglement 
iVe/Coo (Table U). 

The initial simulation cell is a cube of side length Lqj 
where Lq is greater than the typical end-to-end distance 
of the chains. Nch chains are placed in the cell, with peri- 
odic boundary conditions applied in all three directions. 
Each initial chain configuration is a random walk of A^ — 1 
steps with the distribution of bond angles chosen to give 
the targeted large-scale equilibrium chain structure. In 
particular the mean value of cos{6) is adjusted so that 



1+ < cos{9) > 
1- < cos{9) > 



(5) 
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Nch is chosen so that the total number of monomers 
Ntot NNah is 250000 {Lq = 66.5a) for flexible 
{kbend — 0) chains and 70000 {Lq = 43.5a) for semi- 
flexible {khend > 0) chains. The initial monomer number 
density is p = 0.85a^'^. 

After the chains are placed in the cell, we perform 
molecular dynamics (MD) simulations. Newton's equa- 
tions of motion are integrated with the velocity- Verlet 
method [H and timestep St = .007tlj - -OUtlj. The 
system is coupled to a heat bath at temperature T using 
a Langevin thermostat with damping rate 1.0/tlj. 
Only the peculiar velocities are damped. 

We first equilibrate the systems thoroughly at T = 
l.Ouo/fcsj which is well above the glass transition tem- 
perature Tg ~ 0.35uo/kB [1^. The cutoff radius rc 
is set to 2^/^a, as is standard in simulations of melts 
with the bead-spring model [13, H^l • For well-cntanglcd 
chains, the time required for diffusive equilibration is 
prohibitively large. To speed equilibration we use the 
double-bridging-MD hybrid (DBH) algorithm [3l| , where 
Monte Carlo moves that alter the connectivity of chain 
subsections are periodically performed. 

Glassy states are obtained from well-equilibrated melts 
by performing a rapid temperature quench at a cooling 
rate of T = —2 x lO'^uo/ksTLj- We increase rc to its fi- 
nal value, typically 1.5a, and cool at constant density un- 
til the pressure is zero. The quench is then continued at 
zero pressure using a Nose- Hoover barostat 33| with time 
constant IOtlj until the desired T is reached. Larger val- 
ues of Tc lead to higher final densities and larger stresses 
at all strains [1^, but we have checked that using val- 
ues of rc as large as 2.6a does not change the conclusions 
presented below. Indeed, stress-strain curves for different 
rc > 1.5a collapse when normalized by Tfiow Q- In this 
paper T varies from to 0.3uo/fcs. Simulations at T = 
are not directly relevant to experiments, but are useful 
to gain theoretical understanding of polymer deformation 
in the limit where thermal activation is not important. 
To operate in the T — > limit, we remove the Gaussian 
noise term from the standard Langevin thermostat and 
retain the viscous drag term. 

Values of N^, (Table IJ) are measured by performing 
primitive path analyses (PPA) [s^, [131 • Details of the 
PPA procedure are the same as those used for undi- 
luted systems in our recent paper ^ . Melt entanglement 
lengths are consistent with values of Nc from the melt 
plateau moduli ^32,] . Quenching melt states into a glass 
has little effect on the values of Nc determined from PPA 
(SS^. The changes in entanglement density pe = p/^Nc 
upon cooling are primarily due to a 15% increase in p. 
The conclusion that glasses inherit the melt value of iVg 
is consistent with experimental [39j and simulation [30l| 
studies of the craze extension ratio, as discussed in Sec- 
tion llll Al However, it is inconsistent with some entropic 
models of strain hardening that assume that pc in- 
creases rapidly as T decreases. 

The values of N employed in this paper are 12 — 500 
for flexible chains and 4 — 350 for semiflexible chains, 



spanning the range from the unentangled to the fully en- 
tangled {N ^ Nc) limits. It is important to note that un- 
entangled systems {N < Nc) are often brittle. This may 
severely limit the maximum strain that can be studied 
in experiments and complicate comparison to our simu- 
lations. 



TABLE L Chain statistics in fully entangled glasses (N/N^ > 
7) 
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In fundamental studies of strain hardening [3, [3, [i3| , 
compressive rather than tensile deformation is preferred 
because it suppresses strain localization. This allows the 
stress to be measured in uniformly strained systems. The 
rapidity of the quench used here minimizes strain soften- 
ing, which in turn yields ductile, homogeneous deforma- 
tion even at the lowest temperatures and highest strains 
considered. 

We employ two forms of compression; uniaxial and 
plane strain. The stretch Xi along direction i is de- 
fined as Li/L^, where L° is the cell side length at the 
end of the quench. In uniaxial compression, the systems 
are compressed along one direction, z, while maintaining 
zero stresses along the transverse {x,y) directions [4l| . 
In plane strain compression, the systems are also com- 
pressed along the z direction, and zero stress is main- 
tained along the x direction, but the length of the system 
along the y direction is held fixed (Ay = 1). 

Compression is performed at constant true strain rate 
e = Xz/Xz, which is the favored protocol for strain hard- 
ening experiments [40|. The systems are compressed 
to true strains of e = —1.5, corresponding to = 
exp(— 1.5) ~ 0.22, for uniaxial compression and e — —1.2 
(Az ~ 0.30) for plane strain compression. These strains 
are similar to the highest achievable experimentally [Tsj 
in glassy state compression. 

Simulations were performed at e between — 10~^/tz,j 
and —10~^/tlj. As in previous studies of strain hard- 
ening Q and the initial flow stress [3^, \^ , a weak loga- 
rithmic rise in stress with strain rate was observed for 
e| < 3 • lO^"'. This small rise does not change the 
conclusions drawn in the following sections and a sim- 
ilar logarithmic rise is observed in many experiments 
[43| . Thus, while our simulations are performed at much 
higher strain rates than experiments, we expect that they 
capture experimental trends. A more rapid change in 
behavior was observed for |e| > 3 • 10"** and qualitative 
changes in behavior can occur at the much higher rates 
employed in some previous simulations [25l |26| . For ex- 
ample, the time for stress equilibration across the sample 
is of order Lq/cs where is the lowest sound velocity. 
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When this time is larger than the time between plastic 
rearrangements, then each rearrangement occurs before 
the stress field around it has fully equilibrated in response 
to surrounding rearrangements. The decorrelated relax- 
ation of different reg ions leads to a more rapid rise in yield 
stress with rate [35| . The effect of rate on the relaxation 
of unentangled chains is discussed in Section fill Dl 



III. RESULTS 

A. Comparison to Eight-Chain Model 

The eight chain model has been very successful in 
describing the functional form of stress-strain curves and 
is widely used to fit experimental results [28*1. It as- 
sumes that the entanglement network deforms affinely 
at constant volume and employs a body centered cu- 
bic network geometry with eight chains per node. The 
stretch of each segment between nodes is then Xchain = 
+ + A2)/3)1/2, yielding h = \cha^n/ V Ne/C^ in 
Eq. [TJ This choice of network was the main advance 
of the eight chain model over previous entropic models 
[l|, H^l . It allows the model to fit stress-strain curves for 
various strain states, i.e. shear 0, [El a nd uniaxial 0, S] 
or plane strain 0, Q] compression [4^. The prediction 
for the difference between principal stresses along axes i 
and j in the strain hardening regime is 

^.-^.-^;L. + Gi^^^(A?-A^2) , (6) 

where t^i^^ is an independently modeled, rate- and 
temperature-dependent plastic flow stress [H, Q . 

Equation[n] simplifies for the cases of uniaxial and plane 
strain compression considered here and in many exper- 
iments. For uniaxial strain only is nonzero and the 
constant volume constraint implies A^ = Aj, = Aj°'^. For 
plane strain compression, the constant volume constraint 
requires Ax — A^^ and both az and ay are nonzero. 
Equation [5] implies a relation between the strain hard- 
ening of CTz and ay, but the latter does not appear to 
have been measured in experiments. 

Despite its wide use, there are fundamental difficulties 
with the eight chain model that were noted in the Intro- 
duction. As a rubber-elasticity based model, it predicts 
Gr — PeksT. This prediction is about 100 times too 
small at r O.QTg if values of pe are estimated from the 
melt plateau modulus and has the wrong trend with de- 
creasing T d, 0| ■ Some models [1, assume pe is much 
larger than in the melt and rises rapidly as T decreases 
below Tg in order to fit experiments. However, studies of 
crazing in polymer glasses do not indicate any increase 
in Pe over the melt [s^ A constant entanglement 

density is also consistent with the idea that entangle- 
ments represent topological constraints and the observa- 
tion that the topology does not evolve significantly below 
Tg. The extra entanglements added in network models as 



T decreases may capture the effect of glassy constraints 
associated with energy barriers, but it is not clear that 
it is natural to treat such constraints within a rubber- 
elasticity framework. 

Another shortcomin g o f the eight chain model and 
more recent work [4^ |46| | is that the flow stress must 
be introduced as an independent additive constant. Ex- 
periments (47j and our recent simulations Q suggest that 
Tfiow and Gr scale in the same way and are controlled 
by the same physical processes. For example, both de- 
crease nearly linearly as T increases 0, S^j, vanish at 
a strain-rate dependent Tg, and increase logarithmically 
with strain rate (isj . Indeed complete strain hardening 
curves for different rates and cohesion strengths collapsed 
when scaled by tjiow . This suggests that Tfiow is most 
naturally included as a multiplicative rather than addi- 
tive factor. To further test this idea we have examined 
fits to the eight-chain model for a range of fc&end, T and 
strain states. 
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FIG. 1: (Color online) Compressive stress —a^ as a function 
of Az for uniaxial compression at ksT/uo = 0, 0.1, 0.2, and 
0.3 from top to bottom. The chains had kbend ~ l.Siio and 
e = — 10~*/rL,j. Solid lines show a fit to the eight-chain model 
(Eq. inj at fcflT/Mo = 0.2 and the extrapolation of this fit to 
other temperatures using the modified model (Eq. (Tjl. The 
initial elastic rise and yield for 1 > Xz > 0.8 are sensitive to 
aging and are not fit by the eight-chain model. 

Figure [T] shows the compressive stress —a^ as a func- 
tion of Xz for uniaxial compression of systems with 
kbend = I.Smq- Near A^ = 1 there is a sharp elastic 
increase, followed by yield. Both simulations and exper- 
iments [ll m, S Ei, [13] find that this initial region 
(0.8 ^ Az < 1) is sensitive to the past history of the 
sample, including the quench rate and aging. At greater 
compressions the system is "rejuvenated" and the stress 
becomes independent of history . Our discussion will 
focus on this strain hardening regime. 

As with experimental data, the strain hardening re- 
gion (A < 0.8) of all curves can be fit (within random 
stress fluctuations) by adjusting the parameters {rjf^^, 
Gr, Nf.) in Equation [SI The quality of such fits is il- 
lustrated for T = 0.2uo/kB- Typical uncertainties in fit 
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parameters are about 10%. For example, data at all tem- 
peratures can be fit with — 15 ± 1 and best fit values 
lie within this range. Note that the value of — 26 
obtained from the plateau modulus and PPA is substan- 
tially larger 32 1. Fits of Eq. [6] to experiments also 
tend to yield smaller values of Ne than the plateau mod- 
ulus. 

As in previous work 0, M, 113 , fit values of both the 
flow stress and hardening modulus decrease linearly with 
temperature. The temperature where they extrapolate to 
zero, Tg ~ 0.41uo/fc_B, is consistent with previous results 
for the glass transition temperature for this strain rate 
[sst - Data for all temperatures can be fit with a fixed 
ratio P = Gu/Tfiow anywhere in the range from 0.5 to 
0.7. This nearly constant value of /? provides further 
support for the idea that strain hardening scales with 
flow stress. 

The above observations can be incorporated into a 
modified eight chain model that describes the temper- 
ature dependent strain hardening in terms of only four 
parameters t'^io^, Tg, and a temperature independent 
ratio (3. Here tJi^^ is the flow stress at T = and at other 

temperatures Tfiow — 'T^iowi^ ^ '^I'^g)- The shear stress 
in the strain hardening regime is written as 



0,ij 



(7) 



Note that Eq. [7] is equivalent to the usual eight-chain 
model except that it imposes proportionality between the 
flow stress and hardening modulus and assumes a linear 
temperature drop in both. This reduction in parameters 
may be useful in extrapolating experimental data, since 
values for one temperature determine those at any other 
if Tg is known. 

The solid lines in Figure [T] show predictions of Eq. [7] 



based on the fit at T = 0.2uo/kB- t 



flow 



0.634uoa 



(3 = 0.56, and — 14.1. The predictions agree quite 
well with simulation results over an extremely wide range 
ofT/Tg. The largest deviations are of order 10% at T = 
and the smallest Xz- There is a slight over-prediction of 
the change in curvature with increasing T that manifests 
as slight 10%) increases in (3 or decreases in N^, in 
best fits to the data, particularly at T = 0.3. Fits of the 
same quality are obtained for all kbend and strain states 
(see Fig. [5]), suggesting that this simple extrapolation 
may be widely applicable to data from simulations or 
experiment. 

A more stringent test of Eq. [7] is whether it is able to 
predict stresses for multiple strain states with the same 
parameters. As pointed out by Arruda and Boyce 0, 
uniaxial and plane strain compression produce extremely 
different changes in chain configuration. Under plane 
strain compression the chains all stretch in one direction, 
while in uniaxial compression the chains stretch along all 
directions in the plane perpendicular to the compression 
axis. 

Figure [5] shows results for plane strain compression of 
the same systems as Fig. [TJ As before, excellent fits can 



be obtained to Eq. [6] and extrapolations from fits at one 
temperature using Eq. [7] capture the variation in stress 
over the full temperature range. Despite these successes, 
there are troubling inconsistencies in the parameters of 
these fits. While the values of Gr for (Tz in plane strain 
and uniaxial compression are consistent (within our 10% 
uncertainty), the value of is significantly higher for 
plane strain; A^e = 20 ± 2. In addition, the strain hard- 
ening of CTj, and in plane strain are inconsistent. From 
Eq. [SI the value of ay is determined up to an additive 
constant from measurements of Cz. The dashed line in 
Fig. m^b) shows this prediction for T = with the ad- 
ditive constant adjusted to fit data at large A^. At all 
temperatures the variation in (jy is systematically smaller 
than predicted from cr^. The decrease corresponds to a 
reduction in Gr for ay by about 20%. In contrast, best 
fit values of for ay and a^ match to within 1%. 




FIG. 2: (a) Compressive stress — ct^ and (b) perpendicular 
stress — (Jy as a function of Az for plane strain compression 
at ksT/uo — 0, 0.1, 0.2, and 0.3 from top to bottom in each 
panel. The chains had kbend ~ I.Smq and e = — 10~*/rL,j. 
Solid hnes show a fit to the eight-chain model (Eq. [6l at 
ksT /uo = 0.1 and the extrapolation of this fit to other tem- 
peratures using the modified model (Eq. [7| • The dashed line 
in (b) shows the variation in Gy predicted by Eq. [6] and mea- 
sured values of a^. 

A strong dependence of model parameters such as Gn 
and A^e on strain state has also been noted experimentally 
0. One possibility is that the strain state changes the 
role of entanglements in ways that are not captured com- 
pletely by the eight-chain model. Another is that while 
the eight-chain model provides a useful fitting function. 
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it does not capture the correct strain hardening mech- 
anisms. In this case, fit parameters may not have di- 
rect physical significance. The scahng of Gr with a flow 
rather than T supports the view that entropy is not the 
dominant source of stress and this is examined further 
below. It is also important to note that the fits above 
assumed constant volume, but the simulation volume de- 
creased with strain by up to 8% for plane strain and large 
khend- Including the correct Aj in Eq. [6] would change the 
predicted stress by up to 15%. In the following section 
we show that violations of the assumption of affine dis- 
placement of entanglements can produce similar changes. 
Thus the fit parameters compensate for changes that are 
not included in the model, further reducing their direct 
relevance. 



B. Dissipative and Energetic Stresses 

Entropic network models assume that strain hardening 
arises entirely from a reversible increase in the entropy 
of the entanglement network. Experiments showed many 
years ago that strain hardening is also associated with 
large increases in internal energy but this observa- 
tion has not been incorporated into published theoretical 
models. Simulations allow us to separate the role of en- 
tropy and energy in strain hardening. 

For uniaxial and plane strain compression, the stress 
along the compressive axis is directly related to the work 
W done on the system per unit strain: cr^ = dW/ de^ . ctz 
can be separated into an energetic component and a 
thermal component cr^ using the first law of thermody- 
namics: dW — dQ + dU , where U is the internal energy 
of the system and dQ is the heat transfer away from the 
system. This implies 



dW 



dU 
dez 



dQ_ 



(8) 



These quantities are readily obtained from our simulation 
data and could in principle be obtained by differentiat- 
ing results for W and Q from deformation calorimctry 
experiments. Unfortunately existing studies 51, 52, 53] 
have not extended into the strain hardening regime. 

Experimental data are frequently plotted in a manner 
designed to isolate the Gaussian and Langevin contri- 
butions to the strain (Eq. [T|). If is plotted against 
giX) = Xl — Xl, then Eq. [6] predicts a straight fine in 
the Gaussian limit {h << 1). The Langevin correction 
{{3h)~^L~^{h)) adds an upwards curvature. Since A^; is 
not generally measured, experimental stresses are plot- 
ted as a function g{Xz) that is determined by assuming 
constant volume. For uniaxial strain g{\z) — ~ l/Az 
and for plane strain g{Xz) — Xl — 1/A^. 

Figure [3] illustrates the variation of total, thermal 
and energetic stresses with g under uniaxial and plane 
strain compression for the most highly entangled system 
(A^e = 22). The total stress for both strain states shows 
strong upward curvature that is normally attributed to 



the Langevin correction. As expected from entropic net- 
work models^he amount of curvature decreases with in- 
creasing Ne Q. However, Fig. [3] shows that this curva- 
ture is not related to entropy. Almost all of the upward 
curvature is associated with the energetic contribution to 
the stress, while a'^ shows the linear behavior expected 
for Gaussian chains. Similar behavior is observed for all 
Ne and T, and we now discuss the trends in and 
separately. 




FIG. 3: (Color online) Total stress (solid lines) and thermal 
(dashed lines) and potential energy (dot-dashed lines) contri- 
butions for (a) uniaxial compression at ksT/uo = 0.2 and (b) 
plane strain compression at kBT/uo — 0. The systems had 
kt,^„d = 2.0 (Ne = 22), N = 350 and (a) i = -10~^/tlj or 
(b) e = — 10~^/rLj. Dotted hnes show best fits of to Eq. 
□ with (a) Ne = 15.5 and (b) Ne = 22. Straight lines are fits 
to (tJ. Both (Tz and g are negative under compression. 

For all systems, temperatures, and strain states the 
thermal stress is well fit by the linear behavior expected 
for Gaussian chains. There may be a small upwards cur- 
vature, particularly for uniaxial compression, but it is 
comparable to statistical fiuctuations. Attempts to fit 
crj to the eight chain model always require increasing Nf, 
significantly above values obtained from the melt plateau 
modulus. 

We define a thermal hardening modulus Gtherm from 
the slope of finear fits to = ao -\- Gthermg{Xz). Ta- 
ble HI] shows values for Gtherm for uniaxial and plane 
strain compression for various entanglement lengths at 
T = 0.2uo/kB- While Gtherm is systematically higher 
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for uniaxial compression, the differences (10-30%) are 
not large. In contrast, Gtherm decreases rapidly with 
increasing N^- Rubber elasticity theories for Gaussian 
chains would predict Gr cx and it is interesting 

that Gtherm appears to scale in this way. As shown in 
Table HH changes in NeGtherm are within our statistical 
error bars 10%) and show no systematic trend with 

iVe. 



TABLE II: Thermal Moduli Gtherrr^ in units of uo/a^ for T = 
0.2uo/kB and e = — 10~''/rij. Errorbars are about 10%. 
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Figure |4] shows during plane strain compression for 
different N^- Results are shown for T = because the 
energetic stresses are the largest, but results at higher 
temperatures show similar trends. The value of cr^ rises 
to a peak near the yield point, and then drops to a nearly 
constant value for \g\ > 1. The initial behavior for \g\ < 
1 is nearly independent of the entanglement length but 
does depend weakly on age and strain rate. The behavior 
at slightly larger j^l depends strongly on entanglement 
length (Fig. [4|). For example, at T = 0.2uo/fcs the ratio 
of the constant energetic stress to the flow stress rises 
from about 4% for flexible chains to 16% for the most 
entangled system. Similar ratios are obtained for the 
T = data in Fig. H 

Hasan and Boyce examined the enthalpy stored in 
samples of polystyrene (PS), polymethylmethacrylate 
(PMMA) and polycarbonate (PC) as a function of resid- 
ual strain [l^. They found a sharp increase in enthalpy 
up to a strain of about 20% {\g\ = 0.55) that was larger 
in annealed samples than in rapidly quenched samples 
like those used here. For quenched PS, the magnitude of 
the rise in energy density is about 4MPa. Values for the 
work performed are difhcult to extract from the paper, 
but as a rough estimate we take the flow stress (55MPa) 
times the strain (0.2) and find llMPa. Thus in the ini- 
tial stages of deformation, of order a third of the work 
is stored in energy. Calculating the same ratio for our 
simulations gives values between 30 and 45%. 

For strains from -0.2 to -0.8 {\g\ — 2.0) or larger, Hasan 
and Boyce found a weak, nearly linear rise in enthalpy. 
This corresponds to a constant like that observed in 
Fig. |4]for intermediate \g\. Analysis of their figures [l3| 
shows that the ratio of cr^ to the flow stress increases 
from about 4% for PS to 15% for PC. Since PC is more 
entangled than PS, this trend is the same as observed in 
Fig. m Note that in both simulations and experiments 
the fraction of work stored as energy depends strongly 
on the strain amplitude. The fraction stored during the 
initial rise to the flow stress is dependent on sample age 



[T^l and may be of order 50% or more [H, [1^. As \g\ 
increases and remains constant, the fraction stored 
as energy drops towards the much smaller value given by 

Cr^ /cTflow 

Fig. m shows a sharp rise in at the largest values 
of \g\. The onset moves to smaller \g\ and the magni- 
tude of the rise increases as decreases. This rise is 
the source of almost all the curvature in the total stress. 
The experiments of Hasan and Boyce [l^ did not show 
this rise. One reason may be that their experiments only 
extended to uniaxial strains of -0.8 {\g\ = 2.0) for the 
most entangled systems. However, their measurement is 
also limited to the residual enthalpy after the sample is 
unloaded. Simulations were performed to determine the 
amount of energy recovered during unloading from dif- 
ferent strains. The recovered energy is relatively small 
in the region where has a small constant value, but 
rises rapidly at larger strains. Indeed, almost all of the 
energy that contributes to the sharp rise in at large 
\g\ is recovered when samples are unloaded. Thus exper- 
iments could only observe this rise by using deformation 
calorimetry. 




L_, , , 1 , 'j , I , 1 , I , , , L_ 

2 4 6 8 
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FIG. 4: (Color online) Energetic components of stress 
for plane strain compression at T = with strain rate e = 
-W'^/tlj for iVe = 22 (solid line), A^e = 26 (dot-dashed 
line), A^e = 39 (dashed line), and A'^e = 71 (dotted line). 

The rise in energetic stress at large \g\ seems to occur 
when segments of length are pulled taut between en- 
tanglements. To demonstrate this we examined changes 
in chain statistics with stretch. Entropic network models 
assume that the deformation of the entanglement net- 
work is affine to the macroscopic stretch. An affine dis- 
placement would, on average, increase the length of any 
chain segment by a factor of Xchain ■ This stretch can not 
apply at the smallest scales since the length of chemical 
bonds Iq can not increase significantly. As a result chains 
are pulled taut and deform subaffinely on small scales. 
The larger the strain, the larger the length of taut seg- 
ments. 

To illustrate this we compare the rms Euclidean dis- 
tance between monomers separated by n bonds R(n) to 
the affine prediction Raff{n). If i?o(f^) is the distance 
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before stretching, then Raff{n) = XchainRo{n) where 
'^^chain = ('^K + + For the case of plane strain, 

Kha^u = {\l + l + A-2(F/l/o)2)/3, where A. has been 
eliminated by using the ratio V/Vq of the final and initial 
volumes. Volume changes are normally ignored, but are 
large enough to affect the plots shown below. 

Figurein^a) shows the ratio of the observed R{n) to the 
affine prediction as a function of n/Ne for different g{\z) 
and Nf. under plane strain. There is a clear crossover 
from subaffine behavior {R/Raff < 1) to affine behavior 
{R/ Raff — 1) with increasing n. In the subaffine regime 
at small n, chains are pulled nearly taut. The crossover to 
afRne behavior moves to larger n as \g\ increases, imply- 
ing that chains are pulled straight over longer segments. 
For chains with N^^ = 39 the crossover remains slightly 
below iVe at the largest strains considered here. How- 
ever for Ne — 22 the crossover appears to reach by 
\g\ = 5. At larger \g\ the magnitude of R/Raff decreases, 
but the region of rapid crossover appears to remain near 
Ne- This suggests that the entanglements prevent chains 
from stretching taut on longer scales. 

Figure [5{b) shows there is a direct correlation between 
subaffine deformation at N,, and the increase in the en- 
ergetic contribution to the stress. The values of from 
Fig. HI are replotted against R{Ne)/Raff{Ne) instead of 
l^l [5g|. There is a sharp rise in as R{Ne)/ Raf fiNg) 
decreases below about 0.925. As seen in panel (a), this 
corresponds roughly to the point where the length of taut 
segments reaches iVg. This suggests that the energetic 
stress arises when the entanglement network begins to re- 
sist further deformation. As expected from this picture, 
we find a growing tension in covalent bonds as rises. 
However, the maximum tensions in the "worst case" sce- 
nario of plane strain compression for = 22 are only 
about lOOuo/a, which is well below the breaking strength 
240uo/a used in breakable-bond simulations [30j . 

The above findings help to explain some of the discrep- 
ancies in the fit parameters for the eight-chain model. 
The upwards curvature in plots of az vs. \g\ comes from 
energetic terms rather than entropy. Fits to uniaxial and 
plane strain give different Ne because the energetic con- 
tributions are different. However the fit values are never 
far from because the sharp increase in occurs when 
segments of length Ne are pulled taut. Fig. [SJa) also 
indicates that the entanglement network does not de- 
form completely afhnely as assumed in the eight-chain 
model. Even at small strains R{N e) / Raf f {N e) is slightly 
less than one and this would produce significant (^10%) 
changes in the stress from Eqs. [B]or[7| 57]. 

Note that the deviations from affinity observed in Fig. 
[5] are significantly smaller than those predicted from 
rubber-elasticity based models for the non-affinc defor- 
mations in entangled polymers above Tg [58, .59] . These 
models assume that fluctuations about affine deforma- 
tion are confined to the "tube" formed by surround- 
ing entanglements. They predict nonafhne reductions in 
R{Ne) / Raf f (Ne) that are about 50% greater than our re- 
sults at small strains. At the larger strains where we find 



entanglements produce a significant energetic stress, the 
disparity decreases. The discrepancy appears to reflect 
the fact that the nonafflne displacements in our simu- 
lations [§] are much smaller than the tube radius until 
the energetic stress begins to dominate. Interestingly, 
the magnitude of the nonaffine displacements decreases 
slightly with increasing r^, suggesting they are limited 
by cohesive interchain interactions rather than entangle- 
ments. These observations provide further evidence that 
polymers in a glass are not free to explore their tube 
as assumed in entropic models. It would be interesting 
to extend these comparisons to melt models to see what 
additional information can be obtained. 




0.5 1 1.5 2 2.5 

n/N 




R(Ne)/Raff(Ne) 

FIG. 5: (Color online) (a) R{n)/Raff{n) for the same systems 
and conditions as Fig. |4]with n scaled by A^e. Solid lines indi- 
cate A^e = 22 and dashed lines indicate Ne = 39. Curves are 
for \g\ —2.5, 5, 7.5, and 10 from top to bottom, (b) aY plot- 
ted against R/Raff evaluated at n = A^e. This corresponds to 
evaluating a]/ along a vertical slice of (a). Solid, dot-dashed, 
dashed, and dotted lines indicate data for Ne —22, 26, 39, 
and 71 respectively [s^ . 



C. Reversibility and entropic back stresses 

If the work performed in deforming a glass is entropic, 
it should be reversible. However, large strain experiments 
performed well below Tg show only ~10% strain recovery 
upon unloading fl3l| . Entropic network models postulate 
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that there is an entropic "back stress" Q that favors 
further strain reduction but that relaxation is too slow 
to observe because of the high viscosity of the glassy state 
[l|. As expected from this picture, lowering the viscosity 
by heating even slightly above Tg leads to riearly complete 
shape recovery of well-entangled glasses ^ . 

To see if similar behavior occurs in simulations, we 
loaded samples to e = —1.5 at e = —10~^/tlj and T — 
0.2uo/fcs. The samples were then unloaded at the same 
e| and T until all Ci were zero. Fully entangled samples 
(TV = 350, Ne = 39) recovered only 6% of the peak strain 
and unentangled chains {N = 16) recovered slightly less, 
^ 4%. The result for entangled chains is comparable to 
experiments (l3| . 

The samples were then heated to T = QAuo/kB over 
lOOr^j and allowed to relax with a Nose-Hoover barostat 
imposing zero stress in all three directions. Figure [H^a) 
shows the resulting strain recovery. For the entangled 
system, an additional 87% of the strain was recovered af- 
ter IQ^TLj and the rate of recovery remained significant 
at the end of this period. In the unentangled system, 46% 
of the strain is recovered, mainly in the first 2 x 10*tlj. 
While this recovery is substantially smaller than for en- 
tangled systems, network models would predict no strain 
recovery for unentangled chains. Examination of pair 
and bond energies shows that they are nearly constant 
during the relaxation at T = OAuo/kg. These results im- 
ply that entropic stresses drive the relaxation and that 
entanglements play an important role. 

To monitor the entropy in chain confirmations we eval- 
uated the orientational order parameter, P2{cos{a)) = 
(Scos^ia) - l)/2, where cos^{a) =< > / < Rl^ >. 
This quantity measures the deviation from isotropy at 
the end-end scale. There is significant orientation of both 
short and long chains during the initial strain, which is 
discussed further below. During the relaxation to zero 
stress and heating to T = QAuo/kB the orientation re- 
laxes only 3% for entangled chains and 5% for unen- 
tangled chains. As shown in Fig. E^b), rapid and sub- 
stantial deorientation occurs during the strain relaxation 
above Tg. Short chains become nearly isotropic after only 
~ 2 • lO^Tij, and there is little strain recovery after this 
point. Entangled chains deorient more slowly, and both 
P2 and Cz are continuing to evolve slowly at the end of the 
simulations. These results clearly show that the entropy 
of chain orientation drives the strain relaxation. They 
also show that the network of entangled chains prevents 
chains from deorienting without recovery of the macro- 
scopic strain. 

While Fig. [H] provides strong support for an entropic 
back stress, the magnitude of this stress can only be of 
order p^ksTg and thus much smaller than the stresses 
associated with strain hardening. To confirm this we 
took the N — 350 sample studied in Fig. [S] and heated 
to T = 0.4uo/fcs with different stress control. In- 
stead of fixing all Ui to zero, only the total pressure 
p — —{<Tx+(Jy-\-(Jz)/'i was kept at zero while the ratios of 
the Li were fixed at the values after deformation. After 
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FIG. 6: (Color online) Time dependent relaxation at T = 
0.4iio/fcs of (a) true strain and (b) chain orientation param- 
eter P2 for entangled A'^ = 350 (solid lines) and unentangled 
= 16 (dashed lines) systems. Systems were prepared by 
loading to = 1.5 at T = 0.2?io/fcs, unloading to zero stress 
and then heating to 0.4uo/fcs over IOOtlj. 



heating, there was a shear stress — Ox ~ <^z^ f y whose 
direction favored relaxation back to e = 0. The magni- 
tude of this stress relaxed rapidly (~ lO^'r^j) to about 
twice the entropic estimate oi peksTg ~ 0.02uo/a'^, while 
the stress during strain hardening below Tg is more than 
two orders of magnitude larger. Similar results were ob- 
tained at higher temperatures. We next applied a shear 
force of the same magnitude (0.04^0/0"^) to an unstrained 
system at T = 0.4uo/fcs- The magnitude of the strain 
produced by this stress over lO^r^j was 1.2, which is 
comparable to that during stress relaxation (Fig. [5]). 
Both this driven response and the stress relaxation var- 
ied approximately as the logarithm of time, indicating 
that the sample displays creep rather than viscous fiow. 
Recent studies of a similar glassy system also show creep 
behavior at this temperature and time scale [60| . 



D. Chain Length Effects 

The orientation of unentangled chains shown in Fig. 
[SJb) is not expected from entropic network models. For 
N < Nf. there is no entanglement network spanning the 
system. Network models assume that this network is es- 
sential in forcing the deformation of individual chains to 
follow the macroscopic strain. However because chains 
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are not free to relax in the glassy state, chain orientation 
can occur even without entanglements. In recent work on 
uniaxial compression, we found significant strain hard- 
ening of unentangled chains Q and discovered a direct 
connection to chain orientation [l3| as suggested by re- 
cent analytic studies [1^. In this section we extend the 
study of chain length dependence to other strain states 
and systems. 

Figure Et^a) shows stress-strain curves for flexible 
chains (iVe = 71) in plane strain compression for a range 
of N between 12 and 500. At small \g\, the stresses are 
nearly independent of N. Beyond yield, the stresses in- 
crease faster for larger N, reaching an asymptotic limit 
for N ^ Ng as expected from network models Q . How- 
ever, there is significant strain hardening for chains as 
short as iVe/6. Similar behavior is observed for all entan- 
glement densities under both uniaxial and plane strain. 
This is illustrated for /sfcend = 2.0uo {^e — 22) under uni- 
axial strain in Fig. llOf a) and for khend = 0.75mo in Fig. 
2(c) of Ref. [13. 

Examination of individual chain conformations shows 
that strain hardening is directly correlated with increas- 
ing chain orientation [6l|. To quantify this we define a 
microscopic stretch of chains as = Ri/R^ where Ri is 
the rms projection along i of the end-end distance and 
i?° is the value before deformation. Fig. [7] shows for 
flexible chains under plane strain compression. For fully 
entangled chains {N/Ne ^ 7) the microscopic stretch re- 
mains close to the macroscopic stretch as already con- 
cluded from Fig. [5l For A^, the deviation is smaller than 
the line width. For A^ the maximum deviations of about 
2% (lowest dashed line in Fig. EKc)) can be attributed 
to non-afhne deformation of the unentangled ends of the 
chains. 

As the chain length decreases, the microscopic 
stretch shows increasing deviations from the macroscopic 
stretch. Chains compress by less than the imposed strain 
along the z axis, and stretch by a smaller amount along 
the X axis. For each N, A^ is close to the entangled re- 
sults at small \g\ and then saturates at large \g\. The 
onset of saturation in A^ correlates with the saturation 
of the stress, and moves to larger \g\ with increasing N. 
These results clearly show that entanglements force the 
chain orientation to follow the macroscopic stretch but 
that significant chain orientation occurs without entan- 
glements. Strain also orients chains in unentangled melts, 
but is only appreciable when the strain rate is faster than 
chain relaxation times [62]. The extremely slow dynam- 
ics in glasses prevents relaxation of shear-induced orien- 
tation. 

Fig. [ZKd) shows the ratio of the product of the chain 
and macroscopic stretches IliA^/IIiAi. This corresponds 
to the ratio of changes in the volume subtended by the 
chains to changes in the macroscopic volume. For entan- 
gled chains the ratio is close to unity, as expected for a 
crosslinked network. The volume subtended by unentan- 
gled chains need not follow the macroscopic volume, but 
the observed deviations are less than 11% in Fig. \7ld). 



Deviations are even smaller for flexible chains under uni- 
axial strain. 

Figure [8] shows that (Tz is determined directly by the 
microscopic orientation of chains rather than the macro- 
scopic deformation. Results for plane strain compression 
of flexible chains (from Fig. [71Ja)) and uniaxial compres- 
sion of semiflexible chains {N^ = 39) are plotted against 
an effective g calculated from AJ: geff = (A^)^ — (A^)^. 
When plotted against this measure of microscopic chain 
orientation, results for all chain lengths collapse onto a 
universal curve. A similar collapse was obtained in Ref. 
p7j using a single effective orientation parameter X^^^ 
along the compression direction. This was obtained by 
measuring A!^ and using the assumption of constant chain 
volume to determine X^^^ (i.e. X^^^ = ^/^x ^^r plane 
strain). The collapse produced for g{X%f^) is nearly iden- 
tical to that in Fig. [5] because chain volume is nearly 
constant (Fig. [TKd)) and 5e// is mainly determined by 

Xx- 

Fig. [9]shows that results for a^^Oy during plane-strain 
compression also depend only on microscopic chain orien- 
tation. When plotted against the macroscopic \g\, results 
for unentangled chains lie substantially below those for 
entangled chains. When plotted instead against the mi- 
croscopic orientation function ge// = iXz)^ ~ (^^y)^! data 
for all chains collapse onto a universal curve (Fig. [D(b)). 
Note that A^ decreases by as much as 5% from Xy = \ 
for the shortest chains, and this affects the data collapse 

[H. 

The quality of the collapse of the total stress decreases 
slightly as the entanglement length decreases. This is 
illustrated for = 22 {k^end = 2.0) in Fig. [TUl Re- 
sults for fully entangled systems [N > 4iVe) collapse 
completely. Data for smaller N follow the asymptotic 
curve at small |.9e//|, and then drop below it at a |(?e//| 
that decreases with decreasing N. The smallest chains 
in Fig. [TO] and [5Jb) are only a few persistence lengths 
and may not behave like Gaussian chains p^ . However 
such effects are not large enough to explain why results 
for short chains fall below the asymptotic curve in Fig. 
»). 

These discrepancies are instead explained by examin- 
ing the variation with of the energetic contribution 
to the stress. Figure [TOT c) shows plotted against 
g{Xeff). The initial peak at low |5e//| is nearly inde- 
pendent of A^, but the behavior at large |(7e//| is not. 
There is a sharp rise in cr^ for fully entangled chains, 
that does not occur for A^ < 44. The magnitude of this 
rise is comparable to the deviation between results for 
A^ = 44 and the asymptotic curve for entangled chains 
in Fig. [TOTal. These results suggest that while the ther- 
mal contribution to the stress depends only on the chain 
orientation, the energetic contribution at large \g\ only 
occurs for entangled chains. Without the entanglement 
network, chains can contract along their tube to elimi- 
nate the large energetic stresses. 

We have confirmed that increasing the strain rate from 
|e| = 10~^/tlj to does not change the relation 
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between chain orientation and stress described in this 
section. The main effect is to increase Xeff towards A 
with the increase being more pronounced for shorter 
chains. At even higher strain rates, there is almost no 
relaxation, and Xe/f is close to A for all N. This is 
the regime observed in recent simulations [H, HH] with 
atomistic potentials, whose greater complexity requires 
higher strain rates. 



E. Dissipative Stresses and Plasticity 

The large value of the hardening modulus and the 
multiplicative relation between it and flow stress, sug- 
gest that strain hardening is related to dissipation by 
plastic rearrangements rather than stored entropy. To 
quantify the rate of plastic deformation Rp, we examine 
changes in the Lennard-Jones bonds between monomers 
over small strain intervals 6e = 0.005. At the start of 
the interval, all bonds shorter than Tc = 1.5(7 are identi- 
fied. Then the fraction 5f of these bonds whose length 
changes by more than 20% during the interval is evalu- 
ated. This threshold is large enough to exclude changes 
due to elastic deformations. Tests also show that the 
value of Se — 0.005 is small enough that a given atom 
is unlikely to undergo multiple, independent bond rear- 
rangements in any interval. To eliminate activated rear- 
rangements associated with equilibrium aging, the rate of 
plasticity during deformation was monitored at T = 0. 

Figure [TlTal shows the rate of plasticity Rp = Sf/6e 
as a function of |g| during uniaxial compression of fully 
entangled {N — 350) and short {N — 4) chains with 
kbend — 0.75. The two chain lengths lead to very different 
curves, but for both cases Rp is directly proportional to 
the dissipative component of the stress. To illustrate this 
we also plot erf /a* where a* is the constant of propor- 
tionality relating Rp and ct^. Note that even the rapid 
fluctuations with \g\ in the two quantities are correlated. 
These fluctuations are greatly reduced in the total stress, 
which can not be made to correlate as well with Rp. The 
= 4 chains exhibit nearly perfect-plastic behavior for 
\g\ > 1, showing that the correlation is not directly re- 
lated to strain hardening. The fact that a* is nearly the 
same for short chains that flow at a constant Tfiow and 
entangled chains that show significant strain hardening 
at larger strains is clear evidence for the close connection 
between the flow stress and strain hardening. 

In Ref. |17| we showed that Rp and were also cor- 
related for uniaxial compression of chains with ktend — 
and 1.5uo- Figure [T^ shows that this connection extends 
to plane-strain compression. In all cases studied, Rp 
tracks both the mean and local fluctuations. More- 
over, the normalization constants have nearly the same 
value within our numerical uncertainties. Best flts for 
all N, kbend and strain states range between 0.98 and 
l.luo/a^. Since Rp is the rate of rearrangements per LJ 



bond, a* should correspond to the density of LJ bonds 
Plj times the energy dissipated per bond. Each atom 
has on average about 13 L J neighbors and each bond 
is shared by two atoms, so p^j ^ 6.5/3. Thus the en- 
ergy dissipated per bond a* / pLj is about a quarter of 
the binding energy (0.68uo). Note that this value would 
change slightly with the threshold used to define a bond 
rearrangement and other factors in the definition of i?p, 
but the result that a* is similar for all systems with the 
same Tc is more robust. Increasing Vc increases the bind- 
ing energy and also a*. 

The same a* are obtained when the strain rate is re- 
duced to W~^/tlj, but the magnitude of the fluctuations 
increases slightly. Different behavior appears when the 
strain rate is increased to 10~'^/tlj. Fluctuations are 
much smaller since there is insufficient time for stress 
equilibration. There is also a decrease in Rp, while cr^ 
increases. This implies that there are fewer plastic rear- 
rangements involving larger dissipation, presumably be- 
cause the system does not have time to minimize the 
energy. 



IV. SUMMARY AND CONCLUSIONS 

Extensive simulations of strain hardening were per- 
formed for polymer glasses with a wide range of entan- 
glement densities, chain lengths and temperatures. As 
in experiments, we find that the calculated stress-strain 
curves of entangled chains can be fit to expressions de- 
rived from entropic network models (Eq. [S]) . These mod- 
els normally treat the flow stress as an independent pa- 
rameter that is determined by entirely separate mech- 
anisms. However, our simulations and experiments 
47] show that tjiow and Gp are correlated, and that 
both drop linearly to zero as T rises to Tg. This suggests 
that Tfiow enters multiplicatively rather than additively, 
and motivated a simple modification of the eight-chain 
model that describes the full temperature dependence of 
stress-strain curves. While this model may prove useful 
for extrapolating experimental data at one temperature 
to all others, the fit parameters do not appear to have 
physical significance. For example, values of Nf, are dif- 
ferent for uniaxial and plane strain deformation. Another 
difficulty is that the model implies a relation between the 
two nonzero stress components in plane strain compres- 
sion that is not satisfied by the data. We are not aware 
of experimental studies of the transverse stress <jy, but 
it would be interesting to see if the same inconsistency 
could be observed experimentally. 

Separate study of the energetic and thermal compo- 
nents of the stress provided insight into the failures of 
network models. The thermal component of the stress 
scales nearly linearly with \g\ for all systems. Even 
when h is as large as 0.5, the Langevin contribution to 
strain hardening (Eq. [S]) is not evident in cr'^. Instead 
the rapid rise in a at large \g\ is associated with an in- 
crease in the energetic component of stress. This rise 
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is the dominant factor in fits of iVe to the eight-chain 
model. Since the energetic contributions scale in differ- 
ent ways for uniaxial and plane strain, fit values of A^e are 
different for the two strain states. Existing experiments 
have not examined energetic and thermal components of 
the stress separately in the strain hardening regime un- 
der isothermal conditions, but in principle deformation 
calorimetry experiments could do so. Our results sug- 
gest that a study of trends with entanglement density 
would be particularly useful. 

Analysis of chain conformations reveals the origin of 
the energetic stress. As strain increases, chains are 
pulled taut over longer segments. When the length of 
straight segments reaches A^e, entanglements limit fur- 
ther straightening. Additional strain leads to a rapid 
increase in the tension in covalent bonds and the ener- 
getic component of stress. Straightening on a scale of 
order Ng corresponds to large h. Thus even though dif- 
ferent strain states lead to different fit values of Ng in the 
eight-chain model, both fit values tend to follow trends 
in the true Nf,. Modern microscopic theories of rubber 
elasticity (63 | incorporate intra- and inter-chain energetic 
effects due to chain stretching and orientation, respec- 
tively. Analytic studies based on this approach may 
be able to capture the changes in stress and chain con- 
formation observed in our simulations. 

Our simulations reproduce the shape recovery observed 
in experiments when strongly deformed, well-entangled 
glasses are unloaded and heated slightly above Tg [2a |. 
This relaxation is often invoked as evidence for the en- 
tropic stresses predicted by network models. As expected 
from this picture we find a strong correlation between re- 
laxation of strain and the decay of strain-induced orien- 
tational order. However, we show that the stress associ- 
ated with shape recovery is only of order p^ksT and thus 
much too small to account for strain hardening. This 
stress was determined by measuring the shear stress in 
deformed samples after rapid heating, and by identifying 
the shear stress needed to strain an undeformed sample 
at the rate observed in shape recovery. The latter method 
could also be applied in experiments. 

Limited orientation and shape recovery were observed 
for unentangled chains even though entropic models as- 
sume that there is no network to impose chain orienta- 
tion in such systems. Significant strain hardening was 
also found for these unentangled chains. The stress and 
orientation follow results for highly entangled chains at 
small \g\ and saturate at large l^j. The onset of satura- 



tion moves to larger \g\ as N increases. As suggested by 
recent theoretical work [2^ and observed in our recent 
simulations the stress is directly related to effective 
stretches describing the microscopic chain orientation 
rather than the macroscopic stretches . Plots of stress 
against ^(A'^) collapse data for unentangled and highly 
entangled chains onto a single curve. Small deviations 
from this collapse are observed when the energetic con- 
tribution to the stress is large, i.e. when Ng is small and 
strain is large. 

For both entangled and unentangled chains the ther- 
mal contribution to the stress is directly proportional 
to the rate of bond rearrangements. Both the gradual 
trends and rapid fluctuations in the two quantities track 
each other. The proportionality constant is nearly inde- 
pendent of Ne and chain length. The latter result helps 
to explain the connection between Tfio^ and Gr since 
short chains shear at a constant stress near Tfiow There 
has been great interest recently in plasticity in model 
atomic glasses [H, [g^, HI]. It would be interesting to 
check whether the direct correlation between dissipative 
stress and bond breaking/reformation holds in such sys- 
tems [13. 

One of the intriguing questions raised by our results is 
why the thermal contribution to the stress always rises 
nearly linearly with \g\. The stress has the functional 
form expected for the entropy of Gaussian chains even 
in the limit T ^ where entropic contributions to the 
stress must vanish. One possibility is that entropy enters 
indirectly. As \g\ increases, the number of conformations 
available to the chains is reduced. The growing constraint 
on conformations may naturally lead to an increase in the 
number of local bond rearrangements that scales with 
entropy. This would explain the linear increase in the 
rate of plasticity with \g\^ and the corresponding increase 
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FIG. 7: (Color online) Dependence of (a) stress and (b,c) 
microscopic chain orientation on \g\ for flexible chains un- 
der plane strain compression at T = 0.2«o/fc_B with e = 
— 10~^/tlj- Panel (d) shows the ratio of changes in chain 
volume to macroscopic volume. The chains have lengths 
N = 500 (solid lines), N = 107 (dotted lines), iV = 36 (dash- 
dotted lines), A'' = 18 (dashed lines), and TV = 12 (dash- 
dot-dotted lines). The lower dashed line in (c) shows the 
macroscopic stretch A^. 
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FIG. 8: (Color online) Compressive stress as a function 
of microscopic orientation function g^ff = (A^)^ — (X'i)'^ at 
T = 0.2uofcs with e — —10~^/tlj- (a) Plane strain compres- 
sion of flexible chains with lengths A'^ — 500 (solid line), 107 
(dotted line), 36 (dash-dotted line), 18 (dashed line), or 12 
(dot-dash-dot line), (b) Uniaxial compression of semiflexible 

chains (iVe = 39) with TV = 350 (★), 175 (~), 70 (squares), 40 
(•••), 25 (A), 16 ( ) and 10 (o). 
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FIG. 9: (Color online) Stress difference — (jy as a func- 
tion of (a) the macroscopic g and (b) the microscopic orien- 
tation function (?e// = (A^)^ — (A^)^ at T = 0.2uofcs with 
e = -10"Vtlj. Chains have length iV = 500 (solid line), 107 
(dotted line), 36 (dash-dotted line), 18 (dashed line), or 12 
(dot-dash-dot line). 
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FIG. 10: (Color online) (a) Stress as a function of g during 
unicixial compression of kbend = 2.0uo chains with N — 350 
(dotted), 88 (solid), 44 (dash-dot-dotted), 22 (dash-dotted) 
and 11 (dashed) at T = 0.2mo/A:b and e = -10~^/tlj. (b) 
Stress replotted against geff- (c) Energetic stress plotted 
against geff. 
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FIG. 11: (Color online) Rate of plastic deformation Rp (solid 
lines) and normalized dissipative stress (dashed lines) for uni- 
axial compression of ktend ~ 0.75mo chains with A'^ = 350 
(upper curves) and A'^ = 4 (lower curves). Here T = 0, 
e = -IO^^/tlj and cr* = 1.02uo/a^ for iV = 350 and 
l.luo/a^ for iV = 4. 




FIG. 12: (Color online) Rate of plasticity (solid lines) and 
dissipative stress (dashed lines) for plane strain compression 
of flexible (khend = 0) N — 500 chains (lower curves) and 
semiflexible {kbend = 1.5mo) A'^ ~ 350 chains (upper curves). 
Here T — 0, e = —10~^/tlj and a* — 1.05uo/a^ for kbend = 
1.5«o and 0.98uo/q-^ for kbend = Omq- 



